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ABSTRACT 

Context. We consider the radiative transfer problem in a plane-parallel slab of thermal electrons in the presence of an ultra-strong 
magnetic field (B S: B c « 4.4 x 10 13 G). Under these conditions, the magnetic field behaves like a birefringent medium for the propa- 
gating photons, and the electromagnetic radiation is split into two polarization modes, ordinary and extraordinary, that have different 
cross-sections. When the optical depth of the slab is large, the ordinary-mode photons are strongly Comptonized and the photon field 
is dominated by an isotropic component. 

Aims. The radiative transfer problem in strong magnetic fields presents many mathematical issues and analytical or numerical so- 
lutions can be obtained only under some given approximations. We investigate this problem both from the analytical and numerical 
point of view, provide a test of the previous analytical estimates, and extend these results with numerical techniques. 

Methods. We consider here the case of low temperature black-body photons propagating in a sub-relativistic temperature plasma, 
which allows us to deal with a semi-Fokker-Planck approximation of the radiative transfer equation. The problem can then be treated 
with the variable separation method, and we use a numerical technique to find solutions to the eigenvalue problem in the case of a 
singular kernel of the space operator. The singularity of the space kernel is the result of the strong angular dependence of the electron 
cross-section in the presence of a strong magnetic field. 

Results. We provide the numerical solution obtained for eigenvalues and eigenfunctions of the space operator, and the emerging 
Comptonization spectrum of the ordinary-mode photons for any eigenvalue of the space equation and for energies significantly lesser 
than the cyclotron energy, which is on the order of MeV for the intensity of the magnetic field here considered. 

Conclusions. We derived the specific intensity of the ordinary photons, under the approximation of large angle and large optical 
depth. These assumptions allow the equation to be treated using a diffusion-like approximation. 

Key words, acceleration of particles - magnetic fields - radiative transfer - methods: numerical - X-rays: general - stars: magnetars 


1. Introduction 

Regardless of the astrophysical system considered, the presence 
of a magnetic field always plays a fundamental role in defin- 
ing the characteristics and the dynamics of the ongoing phenom- 
ena. Researchers have recently taken steps to include magnetic 
fields in their models of several astrophysical systems, for in- 
stance, in the launching, acceleration, and collimation of rela- 
tivistic jets in gamma-ray burst (e.g. Blandford & Znajek 1977) 
or in the inflow/outflow magnetically-channelled accretion pro- 
cess onto black holes (e.g. Polko et al. 2013) and neutron stars 
(e.g. Thompson & Duncan 1995). 

The increasing interest has dramatically stimulated the de- 
velopment of entire areas of research such as magnetohydrody- 
namics (MHD) to understand how plasma interacts with mag- 
netic fields (see i.e. Lithwick & Goldreich 2001; De Villiers & 
Hawley 2003; Proga et al. 2003). One of the most puzzling top- 
ics related to interaction between magnetic fields and plasma 
is the physics of radiative transfer. Unfortunately an extensive 
and combined treatment of MHD and radiative transfer is still 
an open question, since large numerical simulations are required 
in both fields. 


To reduce the computational time, one possibility is to find 
a satisfactory semi-analytical treatment for radiative transfer 
within a set of reasonable a posteriori assumptions made on the 
basis of observational constraints. Analytical or numerical mod- 
els have the advantage of being much faster than a MonteCarlo 
approach, for instance, and, despite the significantly high degree 
of approximation, provide insights into actual spectra based on 
the proper physical properties. 

What emerges is that the problem is twofold, since the opti- 
cal properties of the plasma are deeply affected by the presence 
of the magnetic field that modifies the interaction between pho- 
tons and electrons through vacuum polarization. In particular, 
when a strong external magnetic field is present, higher order 
quantum effects in the computation of the vacuum polarization 
tensor of a photon within one-loop of a dressed fermion becomes 
important since the strong magnetic field compensates for the 
low value of the coupling constant giving rise to non-linear in- 
teractions among photons, for example fermion-antifermion pair 
creation, photon splitting (Hattori & Itakura 2013). 

In the presence of a plasma with a non-zero density gra- 
dient, the normal modes of photon propagation change from 
being circularly polarized at high electron densities to being 
mostly linearly polarized at low densities. The change in the 
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polarization modes is also accompanied by a change in the 
opacities of the resulting normal modes because vacuum po- 
larization gives rise to significant resonances in the radiation- 
matter interaction. However, for sufficiently high magnetic fields 
( B > 10 14 — 10 15 G) and plasma densities n e k, 10 25 cm -3 and for 
propagating photons of a few keV, both modes are thermalized 
because of the high number of interactions with the plasma. As 
a result, including the off-diagonal vacuum polarization contri- 
butions to the dielectric tensor does not affect the results of the 
transfer calculations (Ozel 2003). 

Within the regime described above, when photons enter a 
magnetized plasma, they split into two polarization modes, de- 
pending on the orientation of their electric field with respect to 

k A B, where k is the momentum of the photon and B is the ex- 
ternal magnetic field. If E || kA B the polarization is called ordi- 
nary; if instead E ±k A B , the polarization is called extraordinary. 

The two polarization modes show different cross-sections 
for the interaction with matter; therefore, a complete solution 
of the problem involves finding the emerging spectra of both 
populations of photons. However, even deriving the exact cross- 
sections for the radiation-matter interaction in the presence of 
strong magnetic fields appeared to be quite demanding, espe- 
cially when the external magnetic field reaches or exceeds the 
critical value B c = 4.413 x 10 13 G so that hv g ~ m e c 2 , where 
v g is the cyclotron frequency which is given by the relation 

eB 

hv R = h = 11.57 Bp keV, (1) 

m e c 

where Bp = B/(10 12 G). The full scattering cross-sections 
were derived in several steps over the years (Ventura 1979; 
Herold 1979; Melrose & Parle 1983; Daugherty & Harding 
1986; Harding & Daugherty 1991). 

Meszaros (1984) discussed the propagation of polarized 
radiation in strong magnetic fields in terms of wave equa- 
tions in the Fourier space considering quantum electrodynamic 
(QED) effects, polarization, and mode ellipticity for Thompson, 
bremsstrahlung, and Compton scattering processes. The depen- 
dence on angles and the frequency of the cross-sections created 
several problems, so the author proposed time-dependent, aver- 
aging, and approximation techniques as a way out of this numer- 
ical impasse. 

Pavlov et al. (1989) presented a comparison of the 
Comptonization process in magnetized and unmagnetized opti- 
cally thick plasmas. Using a simplified analytical treatment and 
including stimulated scattering, they obtained spectra at ener- 
gies lower than the plasma temperature. Later on Lyubarskii 
(1986) and Lyubarskii (1988a, b, hereafter collectively L88), pro- 
posed an extensive analysis of inverse Compton scattering of 
soft photons in a strongly magnetized non-relativistic medium. 
The author considered approximated cross-sections which main- 
tain both the angular and energy dependencies. Using a form 
of the Lokker-Plank approximation which allows the integro- 
differential nature of the radiative transfer equation (RTE) to be 
maintained, he managed to find emerging spectra of ordinary 
photons. 

Lyutikov & Gavriil (2006) proposed a semi-analytical treat- 
ment of the resonant Compton scattering (RCS) process in a 
plane-parallel slab (axisymmetric) configuration. Their model 
consists of the study of the Thomson scattering process of a 
black-body (BB) spectrum in a static plasma filled with an elec- 
tron population at constant temperature and density. Despite 
these crude approximations, this model was successfully applied 
by Rea et al. (2007a, b, 2008) to magnetar spectra below 10 keV. 


Most recently, Nobili et al. (2008) developed a 3D Monte 
Carlo twisted magnetospheric model for the RCS assuming that 
isotropic and unpolarized BB photons are emitted at the neutron 
star surface. 

Nobili et al. (2008, hereafter NTZ08) applied their model to 
a sample of anomalous X-ray pulsars (AXPs) and soft gamma re- 
peaters (SGRs) observed by XMM and INTEGRAL so that the 
energy coverage should be as large as possible. They found that 
the XMM spectra alone are well described by their model, while 
for sources also observed with INTEGRAL (1RXS J 1708-4009, 
1E1841-045, SRG 1900+14) the observed spectra require an ad- 
ditional power-law component with photon index T ~ 1-2. The 
Montecarlo approach used by NTZ08 has the advantage of al- 
lowing the problem to be treated with a sufficiently high degree 
of accuracy, considering both the scattering process (QED ef- 
fects) and the geometry of the problem (multipolar components 
of the magnetic field, non-uniform electron temperature, and 
bulk velocity). 

A detailed analytical and numerical treatment of the radia- 
tive transfer problem in the presence of a relatively strong mag- 
netic field for plasma subjected to inward bulk motion was de- 
veloped by Becker & Wolff (2007, hereafter BW07) who used a 
Fokker-Planck approximation for the RTE (see also Blandford & 
Payne 1981) for the particular case of cylindrical plane-parallel 
geometry. They included a modified angular-integrated form of 
the Thomson cross-section to take into account photon diffusion 
in space and energy over the plasma configuration. However, 
they considered the transfer equation for the zero-moment of 
the photon occupation number, which thus did not take into ac- 
count the possible strong anisotropy in the radiation field specific 
intensity. 

Recently, Farinelli et al. (2012a, hereafter F12) developed a 
numerical code aimed at performing spectral fitting analyses of 
magnetized sources with B <; 10 12 G, like X-ray pulsars and su- 
pergiant fast x-ray transients (SFXT). They adapted a relaxation 
method to seek for solutions of the radiative transfer problem 
in diffusion approximation (see F12 for a detailed discussion of 
the numerical code). The authors successfully apply the model 
to several SFXTs, as reported in Farinelli et al. (2012b). 

Here, we consider the same configuration as L88, i.e. 
a plane-parallel slab of thermal plasma consisting of non- 
relativistic electrons having Thomson optical depth r 0 and tem- 
perature kT e . A uniform magnetic field is oriented along the nor- 
mal of the slab and has a strength of the order of 1 0 14 G or higher. 
A thermal bath of photons distributed as the Planck function with 
temperature A:7bb penetrates the slab from the bottom, acquiring 
a polarization which can be either ordinary (“O” subscript) or ex- 
traordinary (“X” subscript). Vacuum polarization contributions 
and resonant scattering can be neglected considering photon en- 
ergies below 10 keV, magnetic fields of the order of 10 B c or 
higher and electron densities <;10 25 cm -3 . 

In this paper, we reproduce the analytical results obtained in 
L88 and we extend them, introducing numerical techniques to 
solve the RTE of the ordinary photons. 

We compare our numerical results with the case of zero- 
magnetic-field and discuss the main differences between the two 
cases. 

In Sect. 2 we describe the system and the first assumptions 
we made following L88 in order to get a RTE for ordinary pho- 
tons. In Sect. 3 we define the range of angles and optical depths 
within which we can apply the Fourier method to the RTE. In 
Sect. 4 we discuss in details how to solve the equation for the 
space variable, i.e. the optical depth; while in Sect. 5 we show 
how we obtained the solution of the energy equation. In Sect. 6 
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we describe how we calculated the photon angular distribution 
inside the slab and the specific intensity. Section 7 is dedicated 
to the discussion of the results we obtained from the numerical 
code and the comparison with the unmagnetized case. Finally, in 
Sect. 9 we draw our conclusions. 


2. Problem settings 

We consider a plane-parallel slab of non-relativistic electrons 
dipped into a magnetic field B > 4.4 x 10 13 G that is uni- 
form and directed along the normal of the slab. Since charged 
particles follow helicoidal trajectories along the magnetic field 
lines with a gyro-radius which is inversely proportional to the 
magnetic field strength, a one-dimensional Maxwell-Boltzmann 
describes quite accurately the distribution of the electrons. We 
consider non-resonant scattering, therefore photon energies well 
below the cyclotron energy ( hv <K /?v g ), so that the photons have 
two normal polarization modes (Lai & Ho 2002). During the 
scattering process, photons are allowed either to maintain their 
polarization mode or to invert it. The approximated magnetic 
Thomson scattering differential cross-sections for the interaction 
with the plasma are 


dcr 0 ^o 
dcr X ->x (/LaO 

d<To->x (p,p) 
dcr x _,o 



where p = cos i f/ and p! = cos are the cosines of the angles be- 
tween the direction of the magnetic field and the direction of the 
photons before and after the scattering, respectively, while <x T is 
the Thomson cross-section. 

It is worth noticing that only the first term on the right- 
hand side of Eq. (2) does not depend on frequency and on mag- 
netic field through v g , under the assumption of magnetic fields 
stronger than the critical value B c and non-resonant scattering, 
v <s v g . The first term in Eq. (2) dominates over all the other 
ones, except for photons propagating parallel or anti-parallel to 
the magnetic field lines. 

In this regime, the Comptonization process affects the 
emerging spectrum of the ordinary photons, leaving the spec- 
trum of the extraordinary photons almost unchanged. We con- 
centrate then on the study of the radiative transfer process of the 
ordinary photons, although the spectral feature related to the ex- 
traordinary photons can provide an important contribution to the 
total emerging spectrum. The extraordinary photon component 
that will be considered here originates from the fraction of or- 
dinary photons that have changed their polarization during the 
scattering process according to the cross-section (4). 


3. Solution of RTE in a magnetized medium 

The homogeneous integro-differential form of the RTE for 
the ordinary mode photons, neglecting induced processes and 


considering inverse Compton as the leading process, is given by 

(p, v, r) = dp do- 0 ^oi«o (aA v', r) N e (p’) 

- n 0 (p , v, r)N e (p)} - cr 0 ^xN e n 0 (p, v, r), (6) 

where n 0 (p, v, r ) is the ordinary photon occupation number; 
N e (p) is the electron distribution function, which we assume to 
be the one-dimensional Maxwellian; and N e is the electron den- 
sity. For the sake of simplicity, hereafter we will use the simpli- 
fied notations cr 0 o and cr ox for the cross-sections (2) and (4). 

In order to find an approximate solution of the integro- 
differential Eq. (6), we expand both the occupation number and 
the Maxwellian electron distribution in the Taylor series up to 
the second order in Av and Ae, respectively. 

We note that, where an external magnetic field is present, 
even after the Taylor expansion, Eq. (6) does not reduce to a 
purely differential equation both for the space and the energy 
variable (see Eq. (8) in L88), as it does in the general case of 
zero magnetic field (Rybicki & Lightman 1979), hence we are 
still dealing with an integro-differential equation. 

Following the arguments reported in L88, it is worth notic- 
ing that the Comptonization parameter, which in strong mag- 
netic fields is i/nr ~ (2/15)(4k7’ e //;7 e c 2 )T 2 , is smaller than in 
the unmagnetized case by a factor of ~0. 13. Thus, we need to 
assume that our system is optically thick (r » 1), so that the 
average number of scatterings may be large enough to make 
the Comptonization process effective. Owing to the presence of 
the magnetic field, the cross-section (3) approach zero for an- 
gles i // < t -1 , therefore photons travelling at such small angles 
with respect to the field are escaping freely from the plasma 
(Lyubarskii & Sunyaev 1982). 

Clearly, only the photons that move at sufficiently large an- 
gles to the field undergo enough scatterings to be effectively 
Comptonized; in other words, the emerging spectrum will be sig- 
nificantly modified with respect to the seed spectrum. At such 
large angles the optical depth r is large enough to ensure that 
Comptonized photons diffuse almost isotropically. Therefore, 
under the conditions of large angles and large optical depths, 
we neglect the anisotropic part Sn of the occupation number 
n = n i so + Sn ~ «j so . The resulting equation may be handled 
with the separation of variable method. 

We define the dimensionless energy x = hv/kT e and we re- 
place the space variable r with the Thompson optical depth for 
electron scattering r = N e o-jr. We seek a solution of the form 
n(x,r) = s(t)Z(x). Substituting it into Eq. (6) and introduc- 
ing a source distribution accounting for the seed ordinary pho- 
tons A>(x), we obtain the following two independent equations 

f LiZ(x) = A Z(x), 

\ L t s{t) = A s(t ) 


or, more explicitly 

-4 |^;ZM + ZWj - (lx 2 + y) Z{x) = S(x), 

( 1 “ i /l ) = f dr ' K ^ T ~ T 'IW T '), 


where we have defined the quantities 


r = 


15 m t c 2 

TwT' 


and 


/ = 


15 m e c 2 1 
8 kT e x „ 2 


(7) 

( 8 ) 

(9) 
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and the kernel 





(1-p 2 ) 2 - 

e 

P 


Q-" 2 ) 


It-t'I 


dfj. 


( 10 ) 


The quantity A in Eq. (8) is related to the eigenvalue of the space 
operator X T for the eigenfunctions s(r). 

Equation (7) is similar to the standard Kompaneets diffusion 
equation, with no induced processes; however, the term depend- 
ing on the /-parameter on the left-hand side contains the mag- 
netic field dependence, via x g , which is the dimensionless cy- 
clotron energy of the electron, namely x g = hv g /kT e . Although 
the additional magnetic term in Eq. (7) affects the dynamics 
of the Comptonization process, the equation maintains approxi- 
mately the same mathematical structure, so it can be solved with 
the standard Green’s function method (see Sect. 5). 

Equation (8) is instead a homogeneous Fredholm equation of 
the second kind with the logarithmically singular kernel (10). It 
is worth noticing that in the unmagnetized case the space equa- 
tion obtained from the variable separation is described by a dif- 
ferential operator (Rybicki & Lightman 1979). 

The solution of Eq. (8) is not straightforward, since the stan- 
dard integration techniques cannot handle a kernel singularity, 
even if it is moderate. We adopt an algorithm suggested by 
Atkinson & Shampine (2008), described in Sect. 4, which is 
thought specifically for kernels with a quasi-smooth behaviour. 

The overall solution of the radiative transfer problem for 
the isotropic part of the ordinary photon occupation number de- 
scribed by the system (7)— (8) should be found in the form (see 
also Titarchuk et al. 1997; hereafter TMK97) 


n(x, t) = ^ n k (x, t) = ^ c k s k (r)Zk(x), 

k= 1 k=l 


(ID 


where s k (j) is the Mi eigenfunction of Eq. (8) and Zk(x) is the 
solution of Eq. (7) for the Mi eigenvalue. The coefficients c k are 
the Fourier coefficients of the series, i.e. the projections of each 
eigenfunction over a properly chosen spatial photon distribution 

f(r) = e~ T,2T °. 


4. Solution of the spatial problem 

The space problem of RTE (see Eq. (8)) is a homogeneous 
Fredholm integral equation of the second kind, i.e. an eigenvalue 
problem, namely 

L t s(t)= f K(\t - T'\)s(T')dT' = <ts(t), (12) 

J-T0 

where cr = ( 1 - 1/1) and the kernel is given in Eq. (10). Even if the 
kernel has a logarithmic singularity along its diagonal (r = r'), it 
is possible to demonstrate that the integral operator L T maintains 
the property of compactness. Performing the integration over r', 
we obtain an analytical expression for the integrand function 



r'Ddr' 


J>V) 


OV) 


(To-r) 


2-e“ 


dfi. 


(^) 


(to+t) 


(13) 


The integrand is smooth and the integral is finite, so the integral 
operator L T is a compact operator, thus it is bounded and it has 
a complete set of eigenvalues and eigenfunctions (see Atkinson 
1967). 


4. 1 . Numerical treatment of the singularity 

Although, in principle, a logarithmic singularity is integrable, 
we may have numerical problems in mapping the solution in the 
regions of integration near the boundaries. However, the kernel 
can be transformed analytically before the numerical integration 
of Eq. (8) takes place. In particular, defining the variable such 
that t = |r - r'| — » 0 and expanding the kernel around t, we find 



With an additional change of variable p = l/y, it is possible to 
write the kernel as a sum of exponential integrals (Abramowitz 
& Stegun 1964), that for Rett) > 0, are defined as 



The exponential integrals E n (t), if | arg t\ < n, can also be written 
in the form (Bleistein & Handelsman 1986) 


Entt) = 


where 


(-0"" 1 
tn- 1) 


[- log t + <A(«)] - Yj 


n = 0 
m+n— 1 


*A(1) - — i//(n) = -y + 


n~ 1 i 

yi 


(~o m 

(m — n + 1 )m! ’ 


(« > 1 ) 


(16) 


(17) 


and y = 0.57721... is Euler’s constant. The series of exponential 
integrals (16) can be written more concisely as 


Entt)= -P(t) log t+ Qtt), (18) 

where P(t) is a polynomial and Q{t) is a series around t, therefore 
the kernel (10) takes the form 

K(\t - t'|) = ^ [- log(|r - r'|) + Q{\t - r'|)] , (19) 

in which P(|t-t'|) = 1 (Abramowitz & Stegun 1964). In this ex- 
plicit form, the logarithmic singularity has been separated from 
the regular part Q{\t - r'|). 

While the smooth part can be easily treated with any of the 
standard quadrature rules, the integration of the logarithmic term 
requires more attention when t — > 0 (see Fig. 1) where we plot 
separately the behaviour of the two components of the kernel in 
Eq. (19) as a function of the variable t = |r — r'|, with tq = 20 
and the total kernel. A direct integration over the logarithmic 
part of (19) is not straightforward to perform with the standard 
analytical and numerical integration techniques (e.g. Morse & 
Feshbach 1953; Polyanin 2008). 


4.2. The Atkinson-Shampine method 

We adopt the algorithm described by Atkinson & Shampine 
(2008, hereafter AS08). They present a numerical MATLAB pro- 
gram, called Fie, which was created with the purpose of pro- 
viding a numerical code for the integration of Fredholm in- 
tegral equations of the second kind on an interval that can 
be either finite [ a , h] or semi-infinite [0, cx>). The authors con- 
sidered not only kernels K (s, t ) that are smooth functions on 
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t 

Fig. 1. Numerical representation of the kernel (19) with t = |r — r'|. 
Dotted line: singular logarithmic part of the kernel. Dashed line: smooth 
polynomial term of the kernel. Solid line: total kernel, sum of the two 
contributions. 


R = [a, b] x [a, b], but also kernels having a modest singular- 
ity behaviour across the diagonal 5 = t, as long as they conserve 
compactness. 

The integral equations that Fie is called to solve are 
A .r(j) - P log |s - f| x(t ) d? = f(s), a < s <b, (20) 

and they may have solutions x(s) which cannot be necessarily 
smooth at the boundaries of the integration interval. To account 
for the singular behaviour across the diagonal, AS08 uses a prod- 
uct Simpson’s rule for the integration. Hence, AS08 introduced 
a mesh of integration points (to, ■ . . , f„( which is graded near the 
integration limits, a and b, where the behaviour of the solution 
can be critical. The index n is always chosen to be divisible by 4, 
and sufficiently large, in order to guarantee the existence of a 
unique solution of the problem. 

The solution is requested to satisfy the convergence criterion 
for n — » 00 


\\x-x n \U<c\\Kx-K n x\U (c > 0), (21) 

where we have defined the integral operators 

% x = I" log | .v - t\ x(t) dr, (22) 

Ja 

and 'Kn is its approximated form, that we will describe later on. 
Inequality (21) holds if the separation between the mesh points 
is chosen properly. Roughly speaking, the grading of the mesh 
should be intensified near the critical points for the integration. 
In particular, we want the error \\x — x n \\co to be, at least, of order 
0{n~ p ) with p = 3. The general integration scheme, suggested 
by AS08, says that for any triplet of points |r ; _i , tj, tj + \) with j 
odd, the solution x(t) is approximated with a piecewise quadratic 
interpolation function xj(t), so that the integral % x(t) becomes 


% x(t) = I log - t\ jc(r)df = 2_j I logls - t\ x(?)dt 

y=l, 1 
/odd 

n - 1 


r*D ft A r»t 

= I log I*- t\ x(t)dt = ^ I 

d a • 1 J t /_ 

£f 

7=1, J, i - 1 
7 odd 
n 

= ^ M k (s) x(t k ) = x{s). 


log |i - 1 \ Xj{t) dr 


(23) 


where M k (s) are the weights of the interpolating function over 
each subinterval. Thanks to the property of compactness of both 
terms of (19), we use the same mesh to calculate the non-singular 
part of the kernel too and, thus, we define the total weight matrix 
M,(r) = M s k + Mf. 

The solution of Eq. (8) is found by solving the algebraic 
equation that for each eigenvalue <x is the following: 


n 

crx(s) = ^ M k (s)x(t k ) , a < s < b. (24) 

k=0 

It is worth noticing that the actual eigenvalues of (8), A k = 
|(1 - cr k ), are those that appears in the last term of the left-hand 
side of Eq. (7) in the parameter y. 


5. Green’s function of the RTE energy operator 

The energy operator (7) has the form of a confluent hypergeo- 
metric equation and a source term which can be solved with the 
Green’s function method (see e.g. TMK97 and BW07). After 
collecting terms, we obtain a more explicit form of Eq. (7), 
which is 

+ i 4x + * 2 ) % + ( 4x -ix 2 -y)e= s(x x3 x °\ (25) 

where the source term S(x) on the right-hand side has been re- 
placed by a delta function S(x-x o)/x 3 representing a monochro- 
matic source of injected photons. 

Green’s function solution of Eq. (25) can be expressed in 
terms of the Witthaker functions as 


Q (.r, x 0 ) = 


e 2 


»(i- VT+47) 


xq T(2a + 4) 
iTi(a,k, Z,x)e* l( 1+vrs ) 


x 


/ * r +3 

— I(x 0 , a, l), x < x 0 , 

\xol 


i9 r i(a,k,l,xo)e 2 ( I+ ' /1+4/ ) | — j I(x,a,l), x > xq. 


(26) 


where we have defined 

\F\ (a, k, l, x) =1 Fi(a + 2-k,4 + 2a, x Vl +4/) 

\F\ (a, k, l, xq) =1 F\{a + 2 - k,4 + 2a, xq V 1 +4/) 

X °0 2 2 

(.v Vl +4/ + tf + +im f" +1 “ vnS e f dr, (27) 

and the slopes of the two power-laws are determined by the spec- 
tral index 


at = — 
2 


9 15 m e c 2 

— + /IT, 

4 2 kT e 


(28) 


which depends on the kth eigenvalues of Eq. (8). The kth spec- 
trum can be finally expressed as 


f 


Zo.k(x) = Qo (x, xq, A k ) S Uo) dv 0 . 


(29) 
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6. Angular distribution and specific intensity 

The solution for the photon occupation number of the sys- 
tem (7) — (8) for a particular eigenvalue is 

n k (v, t) = s k (T)Zk(y)\ (30) 

thus we can write the specific intensity I = (2/iv 3 /c 2 )n as a se- 
ries of terms that are products of two functions with separated 
dependencies on the variables x and r. Therefore, apart from di- 
mensional factors, the specific intensity is 


I(x,p,r) ~ ^ c k J k (jj,T)Zk(x), 

k= t 


(31) 


where Z(x) is given by (29), in which we dropped the label “O” 
for the sake of clarity. On the other hand, the angular distribu- 
tion is related to the eigenfunctions s k (j) solutions of the space 
problem (8), as described in the relation 


Jkip, t) 


(l V) f e-^- T 'W) 

J-TO 

-(l V)J T °e^ (T '- T W) 


dr' 

A* 

dr' 

A* 


,p > 0 , 

,/r < 0. 


(32) 


The coefficients c k of the expansion of the series (11) are ob- 
tained from the projection of the eigenfunctions over the spatial 
distribution of the source, i.e. 


Clc = f S k (T')f(T')dT, (33) 

J-TO 

where /(r) is a given spatial distribution of the seed photons over 
the bounded medium. 

In addition to the specific intensity carried by ordinary pho- 
tons, since in Eq. (6) we have included also the term which ac- 
counts for mode-switching from O to X, we should calculate the 
contribution to the total specific intensity provided by this popu- 
lation of extraordinary photons. 

The term accounting for the creation of extraordinary pho- 
tons is determined by the cross-section (4), that after the inte- 
gration over i y becomes angle-independent, thus the extraordi- 
nary photons that are originated via mode-switching O — > X, 
can be considered isotropically distributed, with a spectral shape 
(see L88) 

Zx(x) = \ (-) V Zo,k(x) f ° s k (r) dr, (34) 

4 \ X g/ k= i J-TO 

where Zo,k(x) as was calculated in (29). The flux of extraordi- 
nary photons turns out to be quite small with respect to the sum 
of Z k ,o over all k because of the inverse proportionality with 
the magnetic field (x s ~ B ) and the modulation brought by the 
eigenfunctions s k (r). 


7. Results 

In this section we present the results in several steps to discuss 
separately the solutions of Eqs. (8) — (7) and compare them with 
similar results for the same system in absence of magnetic fields. 



Fig. 2. Five series of eigenvalues of Eq. (8) for optical depths r 0 = 
5, 10, 20, 40, 70. Filled dots: 64 eigenvalues for r 0 = 5. Filled squares: 
same set of eigenvalues for r 0 = 10. Filled rhombuses: same set of 
eigenvalues for To = 20. Filled up triangles: To = 40. Filled down trian- 
gles: r 0 = 70. 



Log [ T 0 ] 


Fig. 3. Comparison between the first eigenvalues of the space opera- 
tor Xr defined in Eq. (12) obtained with different methods as a function 
of the optical depth To of the slab. Filled circles: Atkinson-Shampine 
method (Sect. 4). Filled squares: variational approach by L88 (see 
Eq. (36)). Filled rhombuses: asymptotic limit as defined in Eq. (35). 


7. 1. Eigenvalues and eigenfunctions of the space operator 

The algorithm described by AS08, in principle, allows all the 
terms of the infinite series of eigenvalues and eigenvectors of 
Eq. (8) to be found. Nonetheless, the limitation comes from the 
numerical accuracy. 

In Fig. 2 we show five sets of 64 eigenvalues for increas- 
ing maximum optical depth r 0 (from filled circles to filled down 
triangles). The number of eigenvalues that we extrapolate from 
the algorithm is equivalent to the number of points of the grid n 
or less. The limitation that we encounter is numerical and it is 
related to n. 

From Fig. 2 we notice that, above a certain k, the eigenvalues 
start to be indistinguishable, also showing features of numerical 
degeneracy and, as the optical depth increases, the problem mi- 
grates to lower k. Letting n increase, we are able to push the de- 
generacy to higher orders at the expense of computational time. 
With n = 64, we consider eigenvalues up to k = 10 in order 
to avoid this numerical problem. The truncation of the series at 
k = 10, however, does not affect our analysis since, as we will 
show later on in this section, the high energy (Comptonized) part 
of the spectrum is mainly determined by the first term k = 1 , 
while higher orders contribute to the soft peak. 

In Fig. 3, we compare the first eigenvalue that we obtain from 
numerical computations with the analytical estimates performed 
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Fig. 4. First eigenfunctions of the space operator _£ T defined in Eq. (12) 
for different values of maximum optical depth To = 5, 10,20,40,70 
(from solid to dotted line as To increases). 


in L88. Lyubarskii provides two estimates of the first eigenvalue: 
one is obtained by performing a Fourier transform of the ker- 
nel (10), assuming r 0 — > oo, which gives 

2 

d F = ^(log4r 0 -2); (35) 

4r o 

instead, the other estimate, which is 


Avar - 


^(logsr. + r- 



(36) 


is found by solving Eq. (8) with a variational method (here y is 
the Euler’s constant). As suggested by the author, relation (35) 
is no longer satisfied if we are in the case of large optical 
depths (to £ 20). Nevertheless, the eigenvalue obtained through 
Eq. (36) is, by definition, an upper limit of the exact value 
of Ai , hence we expect smaller first eigenvalues for fixed optical 
depths. At small r, the numerical computation and the two ana- 
lytical estimates are very similar, but when the optical depth in- 
creases, the numerical eigenvalues begin to deviate significantly 
from the analytical estimates. 

The left panel in Fig. 4 represents the first eigenfunction for 
different values of maximum optical depth to. The plot shows 
that, in analogy with other physical situations, like a poten- 
tial well with increasing height, eigenfunctions decrease at the 
boundaries, approaching zero as t — > oo. 

Then, we suggest that the larger the optical depth and the 
smaller the number of photons that escape from the slab bound- 
ary, the closer the regime of saturated Comptonization. The re- 
sult is compatible with the analytic expression of the first eigen- 
function considered in L88 (see Eq. (28) of L88II) in which the 
integral term gives a non-zero contribution at t = t 0 . 

In Fig. 5 we compare the spectral index of the emerging 
spectrum for increasing maximum optical depth t 0 for a strongly 
magnetized system in the case of unmagnetized plasma. The 
analytic expression of the spectral index in the former case is 
presented in Eq. (28). Following Sunyaev & Titarchuk (1980, 
hereafter ST80) the spectral index for a slab/disk geometry in 
absence of magnetic field is described by the relation 


OfiJNMAG ~ + A 7 + 


where 


m e c- 


A, 


A = 


12 (to + |) 2 


(37) 


(38) 


Fig. 5. Spectral indexes of the Comptonization spectrum derived from 
the first eigenvalue of the space problem (8) in the presence of a strong 
magnetic field (filled circles) and for the unmagnetized case (filled 
squares). The explicit forms of o-ifro) are given in Eqs. (28) and (37), 
respectively. We note that for a fixed optical depth To of the slab, the 
presence of the magnetic field, through reduction of the scattering cross- 
section, gives rise to softer spectra than an unmagnetized plasma. 


In this case, this eigenvalue also represents the leading term 
in the series (11) and mostly dictates the shape of the emerg- 
ing spectrum. We find that the spectral index of the first 
Comptonization order in the strong magnetic field case is larger 
than in the unmagnetized case for any value of optical depth. The 
magnetic field makes the Comptonization process less efficient 
overall (i/™ ag ~ 0.13 y”°“ as ) as a consequence of the reduction of 
the scattering cross-section (see Eq. (2)) for photons travelling at 
right angles with respect to the magnetic field direction. Unlike 
the case of the radiative transfer problem reported in ST80, the 
presence of the magnetic field induces a significant angular de- 
pendence in the specific intensity, which is taken into account by 
the function J(t,/i) defined in Eq. (32). 

Considering the case of the dominant Comptonization mode 
k = 1 , in the left panel of Fig. 6 we compare the angular dis- 
tribution obtained with the AS08 algorithm described in Sect. 4 
considering the sum of the first ten eigenfunctions and relative 
coefficients and the same distribution, but with the leading term 
only and the angular distribution calculated using the first eigen- 
function corrected at the boundaries (Eq. (29) in L88). 

The distributions appear qualitatively the same, except for 
a scale factor. The normalization gap between the L88 estimate 
and the numerical estimate with the first eigenfunction is due to 
the coefficient c\ multiplying the latter. Basically, the same gap 
is introduced between the two numerical estimates by consid- 
ering the sum of the first ten terms of the series (11) in (32), 
although the shape remains substantially unchanged. Therefore, 
terms with k > 1 contribute to the angular photon distribution as 
a scale factor. 

The right panel in Fig. 6 presents the change in the angular 
distribution for t 0 = 5,20,40. The peak of the distribution be- 
comes narrow and moves towards /j = 1 , where the function goes 
to zero, as the optical depth increases. As expected, for r — > oo 
the function J tends to be flatter, except for /r ~ 1 where the 
peak is located, approaching to an isotropic distribution of the 
photons. 

As for increasing values of optical depth, the probability 
that photons to escape the plasma progressively decreases and 
the eigenfunctions tend to vanish at the boundaries of the slab, 
meaning that only a small number of photons find their way out 
through the plasma. The angular distribution approaches a flat, 
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H 



Fig. 6. Upper panel : angular distribution of ordinary photons for 0 < 
p < 1 and optical depth To = 20. Solid line: angular distribution for 
the sum of the first ten eigenfunctions and relative coefficients c k ob- 
tained with the algorithm described by AS08. Dashed line: same func- 
tion calculated with the first eigenfunction and coefficient from the al- 
gorithm of AS08. Dot-dashed line: same function, but defined as in 
Eq. (29) of L88. Lower panel: angular distribution of the ordinary pho- 
tons emerging at the top of the slab (0 < p < 1) for eigenfunctions and 
relative weights up to k = 10 for the cases of optical depth To = 5 (solid 
line), To = 20 (dashed line), and To = 40 (dot-dashed line). 

isotropic distribution and the system enters the regime of satu- 
rated Comptonization. 



E(keV) 



E( ke V) 


Fig. 7. Upper panel : emergent spectra in arbitrary units, solution of 
Eq. (7) for eigenvalues A k with k = 1,2, 3,4, 5 (see Table 1), maxi- 
mum optical depth t 0 = 20, black-body temperature kT bb = 1 keV, 
and electron temperature kT e = 10 keV. Thin solid line: energy flux for 
k = 1 . Long-dashed line: k = 2. Dot-dashed line: k = 3. Short-dashed 
line: k - 4. Dotted line: k = 5. Thick solid line: total emergent specfi'a 
obtained from the sum of the first five terms of the series weighted with 
the corresponding expansion coefficients c k . The total emergent spec- 
trum has been multiplied by an arbitrary scaling factor in order to be 
compared with its components. Lower panel: emergent spectra in ar- 
bitrary units, solution of Eq. (7), for spectral index a(r, TO defined in 
relation (28) and increasing maximum optical depth with kT bb and kT e 
as specified in the upper panel. Solid line: To = 5. Long-dashed line: 
To = 10. Dot-dashed line: To = 20. Short-dashed line: To = 40. Dotted 
line: To = 70. The values of the spectral indexes and eigenvalues are 
given in Table 2 


7.2. Energy equation and specific intensity 


To find solutions of the energy operator defined in Eq. (7), we 
need to specify a source term on the right-hand side. Here we 
consider the case of black-body seed photons with an exponen- 
tially attenuated spatial distribution, described as 


S ( x , t ) 


Ce - T/2T o yt7 ’3 


X 


3 


s x(kT c lkT bb ) _ 1 ’ 


(39) 


where kTbb and kT e are the photons and the electron tempera- 
tures, respectively. The constant C is a normalization depending 
on the specific problem. 

The energy-dependent part of the specific intensity is ob- 
tained through the convolution of the energy operator’s Green 
function (Eqs. (25)-(29)) with the chosen seed photon distri- 
bution (39). In the left panel of Fig. 7, we show the results for 
the first five term of the series defined in Eq. (29), the solution 
of the problem, i.e. for eigenvalues A k with (k = 1.2, 3,4, 5) 
for the case of to = 20, black-body temperature kTbb = 1 keV, 
and electron temperature kT e = 10 keV. As expected, the high- 
energy part of the spectrum is mostly determined by the first 


Table 1 . First five eigenvalues A k of the space operator (12) and derived 
spectral index a k (see Eq. (28)) for a slab with optical depth To = 20. 


k 

O 

I'- 

ll 

(*k(T = T 0 ) 

1 

0.0085 

0.8460 

2 

0.0297 

2.1922 

3 

0.0555 

3.3498 

4 

0.0838 

4.3623 

5 

0.1130 

5.2496 


Comptonization mode, following the qualitatively similar be- 
haviour of the unmagnetized case (see ST80). 

In Table 1, we present the eigenvalues A k (r = to), found 
as described in the Sect. 7.1, and the spectral index a k for k = 
1, 2, 3, 4, 5 and for r 0 = 20. As we expect from relation (28), the 
larger the eigenvalue, the larger the index, which means that the 
spectrum becomes steeper and steeper. 

In particular, the step between the first and the second eigen- 
value is peculiar; even though eigenvalues and indexes are, in 
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Table 2. Same as Table 1 but for different values of the optical depth r 0 
of the slab. 


To 

^1 (To) 

»l(T 0 ) 

5 

0.0622 

3.6101 

10 

0.0248 

1.9290 

20 

0.0085 

0.8460 

30 

0.0040 

0.4499 

40 

0.0021 

0.2533 

50 

0.0011 

0.1399 

60 

0.0005 

0.0680 

70 

0.0001 

0.0193 


good approximation, equally separated, basically the first term 
provides the most evident deviation from the seed spectrum. 

From a phenomenological point of view, the total spectrum 
can be split into two components: the first term k = 1 of the 
series (11) represents the efficiently Comptonized seed photon, 
while the remaining terms ( k > 2) describe the photons emerg- 
ing from the bounded medium without appreciable modifica- 
tion of their energy. The relative contribution of the two fea- 
tures depends on the spatial distribution of the seed photons (e.g. 
TMK97; Titarchuk & Zannias 1998; Laurent & Titarchuk 1999). 

In the right panel of Fig. 7, we present instead solutions 
relative to mode k = 1 for increasing optical depth r o = 
5, 10, 20, 40, 70 and kT bb = 1 keV and kT e = 10 keV. Varying r 
corresponds to a change in the spectral index a as pointed out in 
Table 2. 

Of course, the larger the optical depth is, the flatter the spec- 
trum becomes. Smaller spectral indexes a imply substantially 
Comptonized spectra for r — > oo up to the asymptotic regime of 
saturated Comptonization. In order to concentrate on the effec- 
tively Comptonized spectrum, we consider the first term of the 
series of Eq. (31) and write the specific intensity for the ordinary 
photons as 

I(x,h,t) ~ ci7i (ji,T 0 )Zi(x), (40) 

where the angular distribution Jiip., to) is calculated in Eq. (32), 
Zi(x) is given by the convolution of Eq. (29) relative to the first 
eigenvalue /li(r), and ci is the first coefficient of the Fourier se- 
ries obtained from the projection of the seed photon space distri- 
bution over the first space eigenfunction. 

We show in Fig. 8 the specific intensity for the ordinary pho- 
tons that keep their polarization during the scattering process, 
by considering the first ten terms of the Fourier series defined in 
Eq. (31) for a given parameter set of seed photons, black-body 
temperature, optical depth and temperature of the electrons of 
the slab, and emerging angle of the radiation field with respect 
to the slab normal. The last variable is mild and does not af- 
fect the spectral shape. Additionally, we also show the ordinary- 
photon specific intensity relative to the first Comptonization 
mode only (k = 1) and the specific intensity of the fraction 
of ordinary photons which become instead extraordinary pho- 
tons. The second feature has been calculated through the rela- 
tion (34) and has been assumed isotropically distributed since 
the cross-section (4) does not depend on the angle of propaga- 
tion of the scattered photons. 

8. Discussion 

The numerical method developed in this paper to investigate the 
spectral formation for photons propagating in a slab dipped into 
a strong magnetic field allows us to consider some qualitative 
and quantitative characteristics of the main spectral features of 
the emerging spectra. 



Fig. 8. Solid line: total specific intensity of ordinary photons obtained 
from the sum of the first ten terms of the Fourier series of Eq. (31). 
Dot-dashed line: specific intensity of ordinary photons for the first 
Comptonization model k = 1 (Eq. (40)). Dashed line: specific intensity 
of extraordinary photons created via mode-switching from a fraction of 
ordinary photons as calculated in Eq. (34). The optical depth of the slab 
is To = 20; the seed photon spectrum is a black-body with kT bb = 1 keV; 
while the electron temperature is kT e = 10 keV. The case of photons es- 
caping with an angle /r = 0.5 with respect to the slab normal is shown. 

Following L88, we have calculated the spectrum of the sub- 
set of the initial seed photon population that has the largest 
probability to be Comptonized, based on cross-sections (2)-(5). 
Considering photons propagating at large angles to the field and 
keeping in mind that we assume photon energies well below 
the cyclotron energy (/zv g ~ 10 MeV for B ~ 10 14 G), ordi- 
nary photons are those which have more chances to undergo a 
considerable amount of scattering. Moreover, the presence of 
cross-section (4) in Eq. (6) tells us that a fraction of the ordi- 
nary photons will change their polarization and will turn into 
extraordinary photons. Separating the variables, we were able to 
solve two independent Eqs. (8) and (7). 

The equation for the space operator (12) is an integral eigen- 
value problem, that we solved numerically by developing the 
technique described in AS07. The eigenfunctions of the space 
operator that we obtain behave like eigenfunctions of many other 
physical systems, like a potential well. As the optical depth 
increases the eigenfunctions approach zero at the boundaries, 
meaning that photons remain trapped in the medium and un- 
dergo a large amount of scattering. The angular photon distri- 
bution calculated from the eigenfunctions of the space Eq. (8) 
reflects the same behaviour. For progressively increasing values 
of the optical depth of the system, the angular distribution flat- 
tens and the photon field becomes almost isotropized. 

The emerging energy spectrum solution of the sys- 
tem (8)— (7) can be essentially split into two main components. 
The first is a soft peak, given by the photons which escaped the 
system without appreciable scatterings, and represented by the 
Comptonization modes ( k > 2) of the expansion series. The sec- 
ond term, corresponding to (k = 1) is the actual Comptonized 
spectrum whose spectral slope depends on the total optical 
depth. For increasing optical depth the spectrum becomes flat- 
ter, approaching the regime of saturated Comptonization when 
t 0 — » oo. Even though the spectral shape and the dependence 
on the eigenvalues is quite similar to the unmagnetized case, 
the comparison between the spectral indexes reveals that for 
fixed physical set-ups, such as electron temperature and optical 
depth, Comptonization is less efficient if a strong magnetic field 
is present. 
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In the calculation of the total spectrum we also took into 
account the contribution of the extraordinary photons created 
via mode switching. However, it provides a small, or even neg- 
ligible, contribution to the overall specific intensity calculated 
here, as we saw in Fig. 8. Nonetheless, the contribution pro- 
vided by the extraordinary photon’s specific intensity is expected 
to be completely Comptonized, since it originates from pho- 
tons that have already scattered several times before the flip in 
polarization. 

It is worth pointing out that the extraordinary mode spec- 
trum obtained in this way does not necessarily represents the 
total contribution of the extraordinary photons to the emerging 
spectrum. It is reasonable to guess, however, that the total ex- 
traordinary spectrum should emerge from the slab basically un- 
changed with respect to their initial distribution, namely a Planck 
spectrum. 

The dominance between the two total contributions, ordinary 
and extraordinary, can be established by several factors, like the 
initial relative percentage of the two seed photon populations and 
the truncation of the series in Eq. (31). Even if we assume that 
both contributions have the same weight in the formation of the 
total spectrum, distinguishing between the two from the point 
of view of data analysis is highly unlikely, especially in the soft 
thermal peak region. 

Therefore, even if solving the coupled system of RTEs 
(Lyubarsky 2002) is the best approach for tackling the problem 
of radiative transfer in strong magnetic fields and describing the 
propagation and the mutual interaction of the two polarization 
modes of the photons, the advantage in the spectral fitting will 
not be so pronounced, since it is very difficult to disentangle the 
two contributions. 

9. Conclusions 

We have performed a numerical study of the radiative transfer 
problem for a plane-parallel slab dipped into a strong magnetic 
field, focusing on the ordinary-mode photons, which have the 
least suppressed scattering cross-section for energies below the 
cyclotron energy, which is of the order of 10 MeV for a magnetic 
field on the order of 10 14 G. The full solution of the RTE is 
obtained by the Fourier method and can be described as a series 
of terms in which the dependencies on the independent variables 
(here energy and optical depth) of the problem are decoupled. 
The simple angular dependence of the Thompson cross-section, 
out of the Klein-Nishina regime, restricts the range of physical 
applicability of the derived spectra up to energies ~100 keV. 

Our work can be considered as a completion and check, 
using numerical techniques, of the analytical results reported 
by L88. A thorough and careful treatment of the singularity of 
the kernel characterizing the space diffusion operator allowed 
us to compute the series of eigenvalues and eigenfunctions of 
the coupled energy-space problem. The implemented numerical 
methods, would in principle give us the possibility of finding as 
many terms as we want of the Fourier series given in Eq. (11), 
even though the ordinary-mode Comptonized spectra are domi- 
nated by the first term. We have compared our solution with the 
analytical estimate given by L88 and we have also performed a 
comparison of the spectral index of the Comptonization spec- 
trum with the case of unmagnetized plasma. 


The geometrical configuration here considered for our com- 
putation, namely a simple static slab with a magnetic field par- 
allel to its normal, makes a straightforward application of the 
model to astrophysical objects difficult, because of the complex 
shape of the magnetic field in magnetar or X-ray pulsars out of 
the neutron star surface (e.g. multipole components). 

In the study of the spectral formation close to the polar caps 
of a neutron star dominated by a dipole component, however, 
our assumptions could work in good approximation allowing us 
to treat the angular dependence of the emerging spectra, at least 
for ordinary-photons, instead of using the zero-moment approx- 
imation for the radiation field. 
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